//
// File: calcAreas.cpp
//
// MATLAB Coder version            : 5.4
// C/C++ source code generated on  : 21-Mar-2025 17:03:22
//

// Include Files
#include "calcAreas.h"
#include "calc_TimeFreq_Ar_types.h"
#include "minOrMax.h"
#include "rt_nonfinite.h"
#include "trapz.h"
#include "coder_array.h"
#include <cmath>

// Function Definitions
//
// calcAreas - Calulates areas/energy under the PSD curve within the freq
// bands defined by VLF, LF, and HF. Returns areas/energies as ms^2,
// percentage, and normalized units. Also returns LF/HF ratio.
//
// Inputs:
//    PSD: PSD vector
//    F: Freq vector
//    VLF, LF, HF: array containing VLF, LF, and HF freq limits
//    flagNormalize: option to normalize PSD to max(PSD)
// Output:
//
// Usage:
//
//
//  Ref: This code is based on the calc_lfhf.m function from Gary Clifford's ECG
//  Toolbox.
//
// Arguments    : const coder::array<double, 2U> &F
//                const coder::array<double, 1U> &PSD
//                const double VLF[2]
//                const double LF[2]
//                const double HF[2]
//                struct2_T *output
// Return Type  : void
//
void calcAreas(const coder::array<double, 2U> &F,
               const coder::array<double, 1U> &PSD, const double VLF[2],
               const double LF[2], const double HF[2], struct2_T *output)
{
  coder::array<double, 1U> b_PSD;
  coder::array<double, 1U> tmppsd;
  coder::array<int, 1U> r10;
  coder::array<int, 1U> r11;
  coder::array<int, 1U> r6;
  coder::array<int, 1U> r7;
  coder::array<int, 1U> r8;
  coder::array<int, 1U> r9;
  coder::array<boolean_T, 2U> r;
  coder::array<boolean_T, 2U> r1;
  coder::array<boolean_T, 2U> r2;
  coder::array<boolean_T, 2U> r3;
  coder::array<boolean_T, 2U> r4;
  coder::array<boolean_T, 2U> r5;
  double aLF;
  double aVLF;
  int b_iindx;
  int c_iindx;
  int end;
  int iindx;
  int loop_ub;
  // normalize PSD if needed
  //  find the indexes corresponding to the VLF, LF, and HF bands
  r.set_size(F.size(0), F.size(1));
  loop_ub = F.size(0) * F.size(1);
  for (end = 0; end < loop_ub; end++) {
    r[end] = (F[end] >= VLF[0]);
  }
  r1.set_size(F.size(0), F.size(1));
  loop_ub = F.size(0) * F.size(1);
  for (end = 0; end < loop_ub; end++) {
    r1[end] = (F[end] <= VLF[1]);
  }
  r2.set_size(F.size(0), F.size(1));
  loop_ub = F.size(0) * F.size(1);
  for (end = 0; end < loop_ub; end++) {
    r2[end] = (F[end] >= LF[0]);
  }
  r3.set_size(F.size(0), F.size(1));
  loop_ub = F.size(0) * F.size(1);
  for (end = 0; end < loop_ub; end++) {
    r3[end] = (F[end] <= LF[1]);
  }
  r4.set_size(F.size(0), F.size(1));
  loop_ub = F.size(0) * F.size(1);
  for (end = 0; end < loop_ub; end++) {
    r4[end] = (F[end] >= HF[0]);
  }
  r5.set_size(F.size(0), F.size(1));
  loop_ub = F.size(0) * F.size(1);
  for (end = 0; end < loop_ub; end++) {
    r5[end] = (F[end] <= HF[1]);
  }
  // Find peaks
  // VLF Peak
  end = r.size(0) * r.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  r6.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      r6[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  end = r.size(0) * r.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  tmppsd.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      tmppsd[loop_ub] = PSD[i];
      loop_ub++;
    }
  }
  //        [pks,ipks] = zipeaks(tmppsd);
  //        if ~isempty(pks)
  //          [~, i]=max(pks);
  //          peakVLF=tmpF(ipks(i));
  //        else
  //          [~, i]=max(tmppsd);
  //          peakVLF=tmpF(i);
  //        end
  coder::internal::maximum(tmppsd, &aVLF, &iindx);
  // LF Peak
  end = r2.size(0) * r2.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  r7.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      r7[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  end = r2.size(0) * r2.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  tmppsd.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      tmppsd[loop_ub] = PSD[i];
      loop_ub++;
    }
  }
  //        [pks,ipks] = zipeaks(tmppsd);
  //        if ~isempty(pks)
  //          [~, i]=max(pks);
  //          peakLF=tmpF(ipks(i));
  //        else
  //          [~, i]=max(tmppsd);
  //          peakLF=tmpF(i);
  //        end
  coder::internal::maximum(tmppsd, &aVLF, &b_iindx);
  // HF Peak
  end = r4.size(0) * r4.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  r8.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      r8[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  end = r4.size(0) * r4.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  tmppsd.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      tmppsd[loop_ub] = PSD[i];
      loop_ub++;
    }
  }
  //        [pks,ipks] = zipeaks(tmppsd);
  //        if ~isempty(pks)
  //          [~, i]=max(pks);
  //          peakHF=tmpF(ipks(i));
  //        else
  //          [~, i]=max(tmppsd);
  //          peakHF=tmpF(i);
  //        end
  coder::internal::maximum(tmppsd, &aVLF, &c_iindx);
  //  calculate raw areas (power under curve), within the freq bands (ms^2)
  end = r.size(0) * r.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      loop_ub++;
    }
  }
  r9.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r[i] && r1[i]) {
      r9[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  tmppsd.set_size(r9.size(0));
  loop_ub = r9.size(0);
  for (end = 0; end < loop_ub; end++) {
    tmppsd[end] = F[r9[end] - 1];
  }
  b_PSD.set_size(r9.size(0));
  loop_ub = r9.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r9[end] - 1];
  }
  aVLF = coder::trapz(tmppsd, b_PSD);
  end = r2.size(0) * r2.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      loop_ub++;
    }
  }
  r10.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r2[i] && r3[i]) {
      r10[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  tmppsd.set_size(r10.size(0));
  loop_ub = r10.size(0);
  for (end = 0; end < loop_ub; end++) {
    tmppsd[end] = F[r10[end] - 1];
  }
  b_PSD.set_size(r10.size(0));
  loop_ub = r10.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r10[end] - 1];
  }
  aLF = coder::trapz(tmppsd, b_PSD);
  end = r4.size(0) * r4.size(1) - 1;
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      loop_ub++;
    }
  }
  r11.set_size(loop_ub);
  loop_ub = 0;
  for (int i{0}; i <= end; i++) {
    if (r4[i] && r5[i]) {
      r11[loop_ub] = i + 1;
      loop_ub++;
    }
  }
  tmppsd.set_size(r11.size(0));
  loop_ub = r11.size(0);
  for (end = 0; end < loop_ub; end++) {
    tmppsd[end] = F[r11[end] - 1];
  }
  b_PSD.set_size(r11.size(0));
  loop_ub = r11.size(0);
  for (end = 0; end < loop_ub; end++) {
    b_PSD[end] = PSD[r11[end] - 1];
  }
  double aHF;
  double aTotal;
  aHF = coder::trapz(tmppsd, b_PSD);
  aTotal = (aVLF + aLF) + aHF;
  // calculate areas relative to the total area (%)
  // calculate normalized areas (relative to HF+LF, n.u.)
  // calculate LF/HF ratio
  // create output structure
  output->aVLF = std::round(aVLF * 100.0) / 100.0;
  //  round
  output->aLF = std::round(aLF * 100.0) / 100.0;
  output->aHF = std::round(aHF * 100.0) / 100.0;
  output->aTotal = std::round(aTotal * 100.0) / 100.0;
  output->pVLF = std::round(aVLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pLF = std::round(aLF / aTotal * 100.0 * 10.0) / 10.0;
  output->pHF = std::round(aHF / aTotal * 100.0 * 10.0) / 10.0;
  aVLF = aLF + aHF;
  output->nLF = std::round(aLF / aVLF * 1000.0) / 1000.0;
  output->nHF = std::round(aHF / aVLF * 1000.0) / 1000.0;
  output->LFHF = std::round(aLF / aHF * 1000.0) / 1000.0;
  output->peakVLF = std::round(F[r6[iindx - 1] - 1] * 100.0) / 100.0;
  output->peakLF = std::round(F[r7[b_iindx - 1] - 1] * 100.0) / 100.0;
  output->peakHF = std::round(F[r8[c_iindx - 1] - 1] * 100.0) / 100.0;
}

//
// File trailer for calcAreas.cpp
//
// [EOF]
//
